Variational calculations for anisotropic solitons in dipolar Bose-Einstein condensates 
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We present variational calculations using a Gaussian trial function to calculate the ground state 
of the Gross-Pitaevskii equation and to describe the dynamics of the quasi-two-dimensional solitons 
in dipolar Bose-Einstein condensates. Furthermore we extend the ansatz to a linear superposition 
of Gaussians improving the results for the ground state to exact agreement with numerical grid 
calculations using imaginary time and split-operator method. We are able to give boundaries for 
the scattering length at which stable solitons may be observed in an experiment. By dynamical 
calculations with coupled Gaussians we are able to describe the rather complex behavior of the 
thermally excited solitons. The discovery of dynamically stabilized solitons indicates the existence 
of such BECs at experimentally accessible temperatures. 
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INTRODUCTION 



Since the prediction and the experimental realization 
of Bose-Einstein condensates (BECs) the field of cold 
atomic gases has been subject of multiple theoretical and 
experimental investigations. BECs with long-ranged in- 
teraction which can experimentally be realized by the 
condensation of atoms with a magnetic dipole moment 
such as ^^Cr are of special interest [IH3]. To stabilize 
such condensates, in general optical traps are applied in 
all three spatial directions. However, it has been shown 
by Tikhonenkov et al. @] that stable quasi-2d solitons 
are possible where a trap is applied only in one direction 
perpendicular to the axis of the aligned atomic dipoles. 

Solitons are a nonlinear effect which arises from the 
dispersion and the nonlinearity cancelling out each other. 
The solitons suffer from two kinds of instabilities: First, 
strong attractive particle interactions can cause the col- 
lapse of the condensate and second, when the interactions 
are only weakly attractive or even repulsive the BEC can 
dissolve in the directions where the external trap is open. 

A prerequisite for the experimental realization of the 
solitons is a detailed theoretical investigation of the pa- 
rameter ranges where the condensate is stable. The rel- 
evant parameters are the trap frequency, the scattering 
length of the contact interaction, and the excitation en- 
ergy which allows for an estimation of the temperature 
range where the solitons can exist. The detailed analy- 
sis of the stationary states and the dynamics of solitons 
based on extended variational calculations is the objec- 
tive of this article. 

In the mean-field approximation the BEC with particle 
number N is described by the extended Gross-Pitaevskii 
equation (GPE). Introducing "natural" units for mass 
TOd = 2m, action h, length ad = (m^oA'^)/(27r?i^), energy 
Ed = /{2ma\), frequency 7d = h/{ma\), abbreviating 
a — flsc/fld and making use of the scaling properties of the 
GPE {f,^,i,i>,E) = {N''^r,N'^-f,N''H,N^/^ti),NE) 
leads to the scaled extended time-dependent GPE in 



"natural" units 
iA^(r) = 



A + 7y + 8^a|V(r)r 



1 - Scos^e 



|^(r')r 



V^(r), (1) 



where the tilde is omitted. The dipole moments of the 
atoms are aligned in z direction by an external magnetic 
field, and 9 is the angle between r — r' and the magnetic 
field axis. In the trap geometry assumed here only a trap 
in the y-direction perpendicular to the magnetic field is 
present, i.e. = Iz = 0. 

In Ref. [1] the GPE ([T]) has been solved approximately 
using a variational approach with a Gaussian type orbital 
and numerically exact by simulations on a grid. The 
grid calculations are numerically quite expensive. For a 
detailed analysis of the parameter space of the quasi-2d 
solitons we therefore introduce and employ an extended 
variational method based on coupled Gaussian functions, 
which has already turned out, for dipolar BECs with an 
axisymmetric three-dimensional trap, to be a full-fledged 
alternative to grid simulations [SHZ]. In this article the 
stable ground state of the condensate is computed by 
imaginary time evolution of an initial wave function, and 
it is shown that typically three to six coupled Gaussians 
are sufficient to obtain fully converged results. For a 
given trap frequency stable solitons exist in a finite 
range of the scattering length a, outside that range the 
condensate collapses or dissolves. 

The dynamics of energetically excited solitons is stud- 
ied by solving the equations of motions for the variational 
parameters in real time. The mean-field energy of the 
ground state is typically only slightly below the energy 
threshold where the soliton can dissolve. However, the 
investigation of the dynamics reveals the existence of dy- 
namically stabilized solitons at energies far above that 
threshold, indicating the possible experimental realiza- 
tion of solitons at temperature T « 5 /xK or even higher 
in the limit of the GPE. The transition temperature of a 
chromium BEC is experimentally given by ~ 700 nK 



2 



[T]. Therefore the soUton can be dynamically stabilized 
in all such BECs. 

The paper is organized as follows: In Sec. |TT] we in- 
vestigate the solitons with the ansatz of a single Gaus- 
sian wave packet. The appealing simplicity of this model 
is that both the stationary states and the dynamics of 
the solitons can be obtained from the Hamiltonian of a 
pseudo particle moving in a three-dimensional potential. 
In Sec. mil the variational ansatz will be extended to a 
linear superposition of Gaussian wave packets (GWPs), 
and the TDVP will be applied to GWPs. The imagi- 



nary time evolution method will be used in Sec. IV for 
the evaluation of the ground state, and the results will be 
compared to the calculations with one Gaussian as well as 
with numerical grid calculations with the split-operator 
method. The analysis of the real time dynamics allows 
us to estimate the stability of excited solitons at finite 
temperatures. Concluding remarks are given in Sec. |V) 



II. VARIATIONAL APPROACH WITH A 
SINGLE GAUSSIAN 

Although the variational ansatz with a single Gaussian 
function cannot provide quantitatively correct results the 
simple model already allows us to gain deep insight into 
the physics of the quasi-2d solitons. 

The ansatz with a single Gaussian as a trial function 
in the TDVP reads 

V;(r) = e<^-'='+^yy"+'^''"+^) , (2) 

with the complex variational parameters 



A, = Al + iAl 



(3a) 
(3b) 



The parameters A^. determine the real half-widths io- = 
1 / \/2A^ of the Gaussian function, 7* describes the am- 
plitude A = exp(— 7') used for normalization of the wave 
function, and ^ is a global phase. The variational ansatz 
(H of course strongly simplifies the problem. Neverthe- 
less the investigation of this model yields physical insight 
and the results for the ground state are, as will be shown, 
qualitatively correct. 

The mean-field energy and the chemical potential are 



-A) 
-A) 



{Vt) 



1 
2 



(4) 
(5) 



respectively. The expectation values with the wave func- 
tion given in Eq. ^ read 



Al 



rx2 



Ai 



r\2 



At 



(6) 




for the kinetic term 
for the trapping potential, 

/ 4* 

for the scattering potential, and 



(7) 



(8) 



(Ki) = \l f ^ {'^.-y Rd (4,4' 1) - 1) ) (9) 



for the dipolar potential with — \J A\l A\. and Hy ~ 



A\IA>- The term 



Rd{x,v,z) = 



dt 



(10) 



denotes an elliptic integral of second kind in "Carlson's 
form" which can be evaluated with a fast and stable ap- 
proximation algorithm [5J [S] more efhciently than the in- 
tegral representation given in [3]. 



A. Hamiltonian form of the equations of motion 

The ground state of the GPE can in principle be ob- 
tained by minimizing the mean-field energy in Eq. 0. 
Alternatively, the TDVP can be applied to derive equa- 
tions of motion for the variational parameters in Eq. (|3| . 
The equations of motion can be used to investigate the 
dynamics of the soliton, and the stationary states of the 
system can be found by searching for the fixed points of 
these equations. The stable fixed point with the lowest 
mean-field energy denotes the ground state. 

For the ansatz of a single Gaussian the system can 
be transformed to the descriptive form of a Hamiltonian 
system by the coordinate transformation 



Pa 



Al 



x,y,z. 



The Hamiltonian in the {q,p) coordinates 



H = 



pI+pI + pI 111 



(11) 



(12) 



IxQyQz 



has the conventional form H = T + Vh- It can be shown 
that Hamilton's equations 



dH 

^= 1^ =P 
op 



P = 



dH 
dq 



dq 



(13) 
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FIG. 1. (Color online) Three-dimensional potential Vh{q) for 
"fy — 2000 visualized by isosurfaces. (a) Scattering length 
a = 0.1. The ellipsoidal form of the isosurfaces marks the 
stable minimum and the hyperbolic form of the isosurfaces 
close to zero marks the saddle point. In (b) the potential 
is rotated to show that the saddle point lies at smaller qy- 
values than the minimum, (c) Same as (a) but at scattering 
length a = 0.08 lowered towards the bifurcation point. The 
minimum and the saddle point approach each other. 



lead to the same equations of motion as the TDVP ap- 
plied to the variational parameters in Eg. ([3| ). 

The potential Vh is visualized in Fig. [l] The isopo- 
tential surfaces close to the stable fixed point have ellip- 
soidal form. At the trap energy i^mf = 7j, the potential 
gets open. The unstable fixed point at smaller values of 



FIG. 2. (Color online) (a) Mean-field energy as function of 
the scattering length for different values of trap frequencies, 
(b) chemical potential and (c) width parameters of the trial 
function for the stable ground state on a logarithmic scale. 



Qy than the local minimum is given by the saddle point 
where the surrounding isosurfaces have hyperbolic form. 
In Fig. [ijc) the potential at a smaller scattering length 
close to the bifurcation point is shown. The local min- 
imum and the saddle point approach one another and 
coincide at the critical scattering length Ocrit- 



B. Stationary states and linear stability 

The fixed points of the system can now be easily cal- 
culated by a nonlinear root search for q = and p = in 
Hamilton's equations (13). In Fig. [2] the mean-field en- 
ergy and the chemical potential of the stable ground state 
and of an excited unstable state are shown as a function 
of the scattering length. The two branches emerge in a 
tangent bifurcation at the critical scattering length Ocrit ■ 
In Fig. the width parameters corresponding to the 
ground state show, that with increasing scattering length 
the soliton gets very broad and finally dissolves. 

To analyze the stability of the solution of the GPE the 
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FIG. 3. (Color online) Real and imaginary part of the eigen- 
values of the Jacobi matrix J as a function of the scattering 
length for the trap frequency 7j, = 2000. The stable branch 
shows pure imaginary eigenvalues, the unstable one has non- 
vanishing real parts. The eigenvalues of the stable As and 
unstable branch A,, coincide at the bifurcation point. The 
inset shows a magnification of the rectangle. 



eigenvalues A of the six-dimensional Jacobi matrix 



d 



x,y, z 



(14) 



A„=0 



or equivalently in the canonical coordinates the eigenval- 
ues — A^ of the three-dimensional real symmetric Hesse 
matrix 9^14/i9g^|vVh=o of the potential Vh(q) are calcu- 
lated. The results are shown in Fig. [3] The eigenvalues 
appear in pairs of different sign and are all pure imagi- 
nary for the stable state. The linear stability for different 
trap frequencies shows qualitatively similar behavior. 



C. Dynamics with a frozen Gaussian 



Hamilton's equations ( 13 ) describe the dynamics of the 



soliton with the potential Vh{q)- The systematic inves- 
tigation and visualization of the dynamics of a Hamilto- 
nian system with three degrees of freedom is a nontrivial 
task. However, the force in the y direction is dominated 
by the strong harmonic trap potential. If therefore 7^ 
is sufficiently large at least for the ground state of the 
condensate Qy takes nearly the value of the harmonic os- 
cillator ground state. As a further simplification of the 
problem we therefore restrict the dynamics to the plane 
given by the condition 



% = Py 



1 

4gf 







1y 



(15) 



which corresponds to a frozen Gaussian ansatz in the qy 
direction. For jy — 2000 the plane qy w 0.011 is marked 
in Fig.[l];b). 
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FIG. 4. (Color online) Periodic trajectories in the potential 
V2D,h at different values of energy for 7^ — 2000 and scattering 
length a — 0.1. The red line denotes the surface of section for 
the PSOSs. There are two periodic trajectories which belong 
to symmetric and antisymmetric oscillations of the soliton. 
The symmetric oscillations are hardly visible lying on top of 
one another and almost parallel to the surface of section. 



The dynamics of the wave function ^ with fixed pa- 



rameters Ay 



0, Al = 1/8q2 



canonical coordinates q^. 
potential 



y — 7j,/2 is described in the 
and qz by the two-dimensional 



V2D,h {qx-.qz 



1 



24V27rg^ 



1 



D 



4:V^qxqz 
2 1 



(16) 



qz ' Mhv ' 



which is illustrated in Fig. |4j 

To describe the dynamics of this two-dimensional sys- 
tem it is convenient to analyze Poincare surfaces of sec- 
tion (PSOS). An adequate choice of the PSOS shown in 
Fig. |4] is using the rotated coordinates and momenta 



cos a sm a 
— sin a cos a 



cos a sm a\ pj. 
— sin a cos a J \pz 



(17a) 
(17b) 



with a — arctan(g2_niin/(j'a;,min) and the crossing condi- 
tion q2 — 0. For constant energy different initial condi- 
tions are integrated. The border of the allowed energy 
range is given by 



(18) 



In Fig. [sja) a PSOS at an energy close to the ground 
state is shown. The motion is completely regular and an 
elliptic fixed point of the Poincare map is present, which 
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FIG. 5. (Color online) Poincare surfaces of sections for the 
dynamics of the two-dimensional potential V2D,h in Eq. (|16|) 



in the rotated coordinates defined in Eq. (171 with trap fre- 
quency 7y = 2000 and (a) a = 0.1, Eyni = 1999.5; (b) 
a = 0.08, E^i = 1995; (c) a = 0.08, E^! = 1999. In 
(a) the motion is completely regular. Lowering the scatting 
length towards the bifurcation point chaotic dynamics ap- 
pears as shown in (b). Increasing the mean- field energy close 
to E^f = 7y at constant scattering length enlarges the chaotic 
regions in the Poincare map, see (c). 



belongs to the antisymmetric periodic oscillation of the 
condensate. The symmetric oscillation is not visible in 
the PSOS for the surface of section being almost parallel 
to it. Increasing the energy towards = jy, the turn- 
ing points of the periodic oscillations in Fig. [4] move to 
larger values of qx and qz (for the symmetric oscillation 
only one turning point, for the antisymmetric both). For 
Emi — ly they lie at infinity, i.e. the soliton dissolves. In 
the Poincare map the resonant tori decay according to 
the Poincare-Birkhoff therorem building the same num- 
ber of elliptic and hyperbolic fixed points in the Poincare 
map. The closer to the bifurcation point the scattering 
length is, the lower the energy is, where this takes place. 



This can be seen especially in Fig. [5] Further increasing 
of the energy leads to regions of chaotic oscillations while 
the elliptic fixed point moves outwards to larger values 
of qi. 

In the picture with a frozen Gaussian the solitons dis- 
solve at energy E^^ — in agreement with [HHO]. If this 
were always true it would imply that in an experiment 
with realistic parameters (e.g., a condensate with 10000 
particles at ^y — 2000 and a = 0.1) the solitons must be 
cooled down to temperatures of about T = 0.15 /zK be- 
cause the energy gap between the ground state and the 
threshold = is v ery small (for a more detailed 



discussion see Sec. IV A). However, a dynamical stabi- 



lization of the solitons at energies E,^[ > jy is possible 
as will be discussed for an ansatz with a single Gaussian 



in the next Section II D and for coupled Gaussians in 
Sec. lrVBl 



D. Three-dimensional dynamics 

In the frozen Gaussian approximation any excitation 
energy of the soliton must be completely deposited in 
the {qx,qz) motion, which leads to the dissolving of the 
soliton at the low threshold Emi — In this section 
we demonstrate that for the fully three-dimensional dy- 
namics with the ansatz ^ a large amount of excitation 
energy can be stored in the qy motion dominated by the 
one-dimensional harmonic trap. 

An example trajectory for trap frequency 7^ = 2000, 
scattering length a — 0.1, and with mean- field energy 
Emi = 3000 far above the threshold at 7y is presented in 
Fig. |6] In (a) and (b) the expectation values 



1 



x,y,z 



(19) 



are drawn as functions of time. The slow oscillations in 
qx{t) and qz{t) generate the Lissajous type motion of the 
quasi-2d soliton visualized in Fig. [6]jc). The fast oscil- 
lation in qy{t) results dominantly from the external har- 
monic trap. The calculations with a single Gaussian thus 
indicate that a rather highly excited soliton may be dy- 
namically stabilized and does not dissolve in contrast to 
the calculations with a frozen Gaussian El [10]. However, 



it must be clarified (in Sec. IV B) whether the dynamical 



stabilization is also possible with coupled Gaussians. 



III. ANSATZ WITH COUPLED GAUSSIANS 

Though the ansatz made in Sec. [H] offers a descrip- 
tive analysis of the soliton, comparison with calculations 
where the GPE is solved numerically on a grid show that 
the results for the ground state only hold qualitatively. 
As it will be shown, the results for the dynamics of the 
soliton are qualitatively different as well. To gain quan- 
titatively correct results the ansatz for the trial function 
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FIG. 6. (Color online) Three-dimensional trajectory of a soli- 
ton at trap frequency 7^ = 2000, scattering length a — 0.1, 
and at mean-field energy E^f = 3000 far above the threshold 
Emf = 7j/ where dissolving is possible, (a) and (b): Time de- 
pendence of the expectation values qcr{t) with a — x,y, z. (c) 
Lissajous type motion of the projection in the {qx,qz) plane. 



has to be modified. As shown in |6l [7] the condensate 
wavefunction can be well described by a superposition 
of Gaussians. We therefore apply the TDVP to coupled 
Gaussian wave packets [TTI - IT^ . The ground state will be 
determined by the imaginary time evolution (ITE) with 
the norm conservation embedded in the TDVP as a con- 
straint. 

In this section we elaborate the theory and derive the 
formulae for the computations with coupled Gaussians. 
The results for the stationary states and the dynamics of 
the solitons are presented in Sec. llV) 



A. Time-dependent variational principle for GWPs 

In this article the TDVP provides the basis for the 
solution of the time-dependent GPE. Its application to 



GWPs was originally introduced by Heller [TSIIIS] for the 
description of atomic and molecular quantum dynamics 
and later on applied to BEG [51 [7] . For the convenience 
of the reader we here shortly recapitulate this method. 
The quantity 



/ = - Hx{t)\\^ = min 



(20) 



is to be minimized with respect to cf). Afterwards it is 
set (j) = X- The wave function shall be parameterized 
by the variational parameters z(t). The variation of / 
in Eq. (20) carries over to the variation of z and z and 



results in 



SI 



{Sx\ X) 

dx 
dz 



Sz 



(Xl Sx) - mi Hx) - {Hx\ Wx> 
dx 



X + '-^Hx 



X + ^Hx 



dz 



Si 



(21) 



for the wave function x being constant. Because the vari- 
ational parameters are complex quantities both brackets 



in Eq. (21 ) have to vanish separately. This yields 

Hx 



dx 
dz 



oz 



dx 
dz 



which can be written shortly as 

K ■z = -ih. 



(22) 



(23) 



where K is an hermitian positive definite matrix. The 



linear system ( 23 1 has to be solved for every time step to 



integrate the equations of motion i — J{z). 

We now choose a linear superposition of N Gaussians 
as the trial function 



N N 



k=l 



k=l 



N 

fc=i 



where A'' — A'''^ + iA'^'^ are complex diagonal matrices of 
dimension 3x3 and 7'^ = 7'=''' -|-i7'='* denotes the relative 
phases and the amplitude of the Gaussians, respectively. 
The Gaussians are fixed at the origin. The time evolution 
can be considered as the motion in an effective time- 
dependent harmonic potential 



VcS (x) = Wo + ^XV2X . 



(25) 



The dynamics of the GWPs are now determined by 
the TDVP, which variationally fits the effective time- 
dependent harmonic potential coefficients Vq,V2 to the 
underlying potential. Splitting the Hamiltonian H = 
T+V and operating on the trial wave function (24) yields 



N 

i^ _ = ^ + 2i Tr A'' 



k=l 



+ x(-A^ ^A{A^f)x\g{z^,x). (26) 
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The parameters Vq are scalars and the matrices V2 are 
diagonal matrices. The equations of motion then read 



7'= = 2iTrA'= - , 



(27a) 
(27b) 



has to be calculated which yields the system of equations 

(34) 



K 




M 







with 



The coefficients Wg and V2 have to be calculated from 
the TDVP 



K = 



i* - 




h = 




(28) 



E f^o + \xV^x\ g{z\ x))=0. (29) 



k=l 



In the special case of Gaussian wave packets the equa- 
tions of motion (split into real and imaginary part 
and combining the variational parameters to z — 

be written in the form 



Combining the coefficients Vq and V2 in a complex vector 
V the set of equations ( 29 ) can be written as 



z ^Uv + d, 



(35) 



Kv 



(30) 



The positive definite hermitian matrix K includes in con- 
trast to Eq. (23) the kinetic operator. Eq. (30 1 has to be 



where U consists of the terms linear in v and the con- 
stant terms are absorbed in the vector d. The constraints 
combined in / = (/i, . . . , /„) yield 



solved for every time step and the resulting vector v has 
then to be inserted in Eq. (27 1. 



oz 



5/ , 
oz 



Uv + d^O. 



In the numerical integration of the equations of motion The set of linear equations (|34|) then reads 



(271 (e.g. with a standard algorithm like Runge-Kutta) , 
the quadratic term can lead to the numerical overflow 
[T5] . It is therefore reasonable to split up the matrices A'^ 
into two matrices C and B according to A'' — [C^)^^ . 
For the case of diagonal matrices A^ , and B^ are di- 
agonal, too. Hence the splitting can be done component- 
by-component. This leads to 



B'' = -^C'^V^ , 
2 ^ 



(31) 



where C(0) = 1 and B(Q) = A{0) are taken as initial 
values. 



(36) 



(37) 



(38) 



Solving this system of linear equations and inserting the 
result vector in Eq. (27) enables one to integrate the 
equations of motion. 



with 




B. Constraints in tlie TDVP 



C. Energy functional 



The most important reason for the introduction of con- 
straints is the conservation of the norm in the imaginary 



time evolution (cf. Sec. HID I. Apart from that inequaf 



ity constraints can be used to avoid matrix singularities 
which arise from overcrowding the basis set |17| . In prin- 
ciple m constraints can be introduced combined in the 
771-dimensional vector /. The constraints are embedded 
then by a set of Lagrangian multiplicators A € R™ 



L = 1 + XMz . 



(32) 



where M = §C is a real m x 2n„ matrix. To find the 



dz 

minimum of / 



dL 



with a; 



(33) 



To gain the mean-field energy of the soliton the time- 
dependent variational parameters of the GWPs are used 
in the evaluation of the integrals in the energy functional. 
The energy functional for N coupled Gaussians reads 



1 

+ 2 



Lk 



E 

Lk 



Ai* - AS 



(39) 



With (7 = x,y,z the integrals in (39) yield 
(*|-A|M/)=E(.9'hA|/) 



(40) 
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(*|14|*) = -JiE 



(.91 5'), (41) 



n 



Ak I Ai — Al* ^ A^' 



,i(7'''+7'-7'*-7^*) 



(42) 



i.j,k.l 



n 



with the abbreviations 



^ — ^ (ylj — A!^ ) (^A], + — yli — Ai, ^ 



and 











{A^-Ai') 




Al — Ai 


{^y ^ ^y ) 




A^: -Ai*) 



The additional integrals in the TDVP {a'^V{x)), 
(a'^V{x)) and l^a^ fi'^V [x)) are obtained easily from the 
integrals in the energy functional. 

D. Imaginary time evolution 
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FIG. 7. (Color online) Imaginary time evolution of 10 coupled 
Gaussians. The parameters used are 7^, = 2000, a = 0.1. 
Plotted are, from the upper to the lower panel, the real part 
of the width parameters, the imaginary part of the width 
parameters, real and imaginary part of 7 and the mean-field 
energy. The mean-field energy converges to the ground state 
of the system as described in the text. For comparison the 
mean-field energy in a fast calculation with coarse integration 
accuracy is shown as well. 



In principle the ground state of the system of coupled 
Gaussians could be obtained by a nonlinear root search in 
the equations of motion. Since AN complex initial values 
are required, it is difficult to find the fixed points this way 
for a large number of Gaussians N . The imaginary time 
evolution represents an alternative for the determination 
of the ground state. 

In the linear Schrodinger equation the transformation 
t ^ T = \t leads to an exponential damping of the excited 
states and converges to the ground state for r — > 00. The 
method can be applied to the nonlinear GPE as well. In 
imaginary time the norm is not conserved any more. This 
is crucial in the nonlinear GPE because the decay of the 
norm has the meaning of particle losses. To avoid the 
norm decay we implement the norm conservation as a 
constraint 



k,l 



(44) 



fkl 



in the TDVP as elaborated above. Supposing a general 
complex rotation of the time r = e'^t with ^ = -|- 



i^' = — e the equations of motion can be written in 
the general form 



z = Uv + d 



with 



U 





(0)3^> 




)nxN 


i^)3NxN 






<3N 


(^)nx3N 


)3Nx3N 


)nxN 






i^^)NxN 


i^)3NxN 


\i^)Nx3N 


(e)3^> 


<3N 


(0)jVx3Af 


ii^)3Nx3N 



(45) 



2 Tryl'^^* 
-2 TtA''''' 

2 Tr A''-'' 

2 Tr A'''^ 
\A {A'^^'^f -4{A''^'yj 



V = 



9 



where the abreviative notation in U can be understood 
as foUows: Every entry (A)axfc is an a x 6 submatrix [/'. 
If A = all elements of U' 0, otherwise C7' = A • 1. The 



vectors Vn — 





l,r 



I It 

rN 



N,r\'j 



, V, 



N,i\T 



N,i\T 
) 



con- 
tain the real and imaginary parts of the coefficients of the 
effective potential. The entries in d are vectors with the 
index k running from 1 to fc. One special case is the 
real time evolution = 1, = and the other special 
case the imaginary time evolution = 0, = 1. The 
gradient vector — df /dz in Eq. ( 36 ) is given by 



df 
df 
df 
df 



N 



2^Im/fcz, 

1=1 

N 

^Im 

1=1 "''^ 

N 

2^Re/, 



fkl 
,kl 



1=1 
N 



(46a) 
(46b) 
(46c) 
(46d) 



1=1 



with which one obtains U — ■ U and d — F'^ ■ d. 

With this method the norm is conserved during the 
numerical integration for rather long times. However for 
very long times a small drift in the norm is present due 
to the numerical error of the integration. The ITE of a 
system of 10 coupled Gaussians is shown in Fig. [7| 

If the relative accuracy of the integrater is cho- 
sen coarse, the mean-field energy does not decay 
monotonously but for long times it converges to the same 
value as an integration with a fine relative accuracy thus 
a coarse accuracy can be chosen for the sake of time. The 
ITE is very robust to the initial choice of the wave func- 
tion. This also holds for large numbers of Gaussians cou- 
pled, where the nonlinear root search often fails. How- 
ever, only the stable ground state is accessible with this 
method. 
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FIG. 8. (Color online) (a) Mean-field energy as a function 
of the scattering length of the ground state of a soliton at 
trap frequency 7^ = 2000 obtained by the ITE with different 
numbers of Gaussians. For comparison the result for a single 
Gaussian, obtained by a nonlinear root search and the results 
obtained by numerical grid calculations are plotted as well. 

(b) The chemical potential for the same parameters as above. 

(c) Temperature for different particle numbers belonging to 
the difference of upper limit E-^i — "fy and ground state en- 
ergy 



A. Stationary ground state 



IV. RESULTS WITH COUPLED GAUSSIANS 

Applying the variational ansatz with coupled Gaus- 
sians to the wave functions of dipolar condensates we are 
now able to obtain quantitatively correct ranges of the 
parameters where stable quasi-2d solitons can exist. For 
an experimental realization of solitons with chromium 
atoms realistic trap frequencies are typically a few hun- 
dred Hertz for condensates with about 10000 to 20000 
particles. These values roughly correspond to an interval 
of 2000 < 7y < 20000. We therefore focus our calcu- 
lations on two values for the scaled trap frequency, viz. 
7y = 2000 and = 20000. 



In Figs. |8] and [9] the results for the ground state of 
the soliton are shown at trap frequencies 7^ — 2000 and 
7.y = 20000, respectively. In Figs. [8][a) and [9]ja) the 
mean-field energy as a function of the scattering length 
is plotted for up to 6 coupled Gaussians. For compari- 
son the calculation using a nonlinear root search for one 
Gaussian as discussed in Sec. |IIB| is displayed, as well. 
The mean-field energy and also the chemical potential 
[see Figs, [sjjb) and [9]jb)] converge with increasing num- 
ber of Gaussians N. The detailed convergence proper- 
ties of the calculations with increasing number of Gaus- 



sians are illustrated in Fig. 10 where calculations for up 
to 12 Gaussians are shown at constant scattering length 
a = 0.1. The ground state energy converges fast with the 
number of Gaussians and the corrections above a number 



10 



20000 
19800 
19600 
19400 
19200 
19000 



. ' ' ' ' ' .... ' 




^<^' ' (a) - 




root search 1 Gaussian 




ITE 1 Gaussian 




ITE 2 Gaussian 




ITE 3 Gaussian 




ITE 6 Gaussian 




numerical srid calculation X 




, £;f=Yy = 20000^ 




a 6 
H 

4 

2 






10000 particles -- 






20000 particles 

















FIG. 9. (Color online) Same as Fig. [s] but for trap frequency 
7j, = 20000. Compared to the results for 71, = 2000 the 
range of the scattering length is smaller but the temperature 
is higher. 



of about 5 Gaussians can be neglected. The small energy 
differences between the variational and the grid calcula- 
tions in Fig.|9]^a) may be due to the finite grid size, which 
is limited by an acceptable computation time. 

In the calculations with a single Gaussian a stable and 
an unstable state are created in a tangent bifurcation at 
a scattering length a = Ocr- Using coupled Gaussians the 
critical scattering length where the condensate collapses 
is shifted to higher values, e.g., from a^=^ = 0.0734 to 
a^"^ = 0.0820 at trap frequency — 2000 and from 
a^=i = 0.1197 to a^=^ = 0.1246 at trap frequency 

— 20000. The ITE can only provide the stable ground 
state but no unstable states which are certainly also in- 
volved in the bifurcation. For dipolar condensates in an 
axisymmetric trap a complicated bifurcation scenario has 
been revealed [U [7] , and it will be an interesting future 
task to study bifurcations of the soliton states in more 
detail. 

The corrections to the calculations with a single Gaus- 
sian decrease with growing scattering lengths. This can 
be understood when looking at the spatial size of the 
soliton [cf. Fig. 2]jc)]. For large expansion of the con- 
densate the nonlinear interaction terms in the GPE are 
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FIG. 10. (Color online) Mean-field energy for trap frequency 
yy = 2000 and fixed scattering length a = 0.1. The mean-field 
energy converges fast with increasing number of Gaussians. 



small and thus the mean-field energy is dominated by 
the ground state energy of the harmonic trap in the y 
direction. Above a certain threshold value for the scat- 
tering length the soliton is no longer bound but dissolves. 
The thresholds at, e.g., a = 0.137 for jy — 2000 and 
a = 0.157 for 7y = 20000 obtained with a single Gaus- 
sion (see Fig. I2| arc only very slightly shifted to higher 
values when using coupled Gaussians. 

The extended variational approach allows us to com- 
pute for each trap frequency an accurate lower and upper 
critical value of the scattering length, i.e., a range where 
stable solitons in dipolar BECs can exist. However, in an 
experiment thermal excitations of the soliton may cause 
the dissolving of the soliton, and therefore it is also neces- 
sary to determine an energy or temperature limit for the 
existence of stable solitons. As a dissolving condensate 
without any interaction must have at least the zero point 
energy = jy of the harmonic trap in the y direction 
the difference between this threshold and the mean-field 
energy of the ground state. 



Iv 



E: 



mf 



mkeT 
2E. 



(47) 



can be taken for a very conservative estimation of the 
temperature at which the soliton should survive ther- 
mal excitations. In Figs.jsjjc) and[9]Jc) that temperature 
is plotted as a function of the scattering length for two 
condensates with N = 10000 and N = 20000 particles. 
With this conservative estimate a soliton at, e.g., trap 
frequency jy = 2000 and scattering length a = 0.1 must 
be cooled down to about T — 0.15 ^K. However, soli- 
tons may exist at much higher temperatures when they 
are dynamically stabilized as discussed in Sec. |II D| for 
the ansatz with a single Gaussian. We now extend that 
discussion to the ansatz with coupled Gaussians. 
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B. Dynamics with coupled Gaussians 

The ansatz with a single Gaussian used to describe the 
dynamics of the sohtons in Sec. |IID| and especially the 
further simplification with a frozen Gaussian in Sec. |II C| 
are basically mathematical model systems, and the re- 
sults obtained with these models have to be seen from an 
academic point of view. For a realistic description it is 
necessary to investigate the dynamics of the solitons with 
the extended ansatz ( 24 1 of coupled Gaussians which has 



shown to significantly improve the results for the ground 
state in Sec. |IV A[ Of special interest is to search for the 
existence of dynamically stabilized solitons at energies 
above the threshold E^f — To this aim the equations 
of motion ( 27 1 are integrated with the initial wave func- 



tion chosen appropriately to yield the desired mean-field 
energy. 



The condensate wave function (24) with N coupled 



Gaussians is parametrized by 8iV real time-dependent 
variational parameters. For the visualization of the wave 
function we use the reduced set of parameters 



(7 = x,y,z 



(48) 



which are related to the extension of the soliton in the 
three spatial dimensions and can be directly compared 
with the results for a single Gaussian using the coordi- 
nates in Eq. (19). 



For energies close to the ground state energy peri- 
odic oscillations can be found. The parameters of the 
turning points of the periodic oscillations show a rather 
complicated behavior when the mean-field energy is in- 
creased. The different periodic oscillations vanish in mul- 
tiple bifurcations close to the ground state energy. Above 
^-mf = ly no periodic oscillations were found. However 
increasing the energy above i^mf = ly yields oscillations 
of the soliton which do not destroy it. Applying a frozen 
Gaussian approximation yields stable oscillations only 
slightly above this limit. But allowing oscillations in all 
spatial directions stabilizes the condensate and gives rise 



to an oscillating soliton as shown in Fig. 11 for a tra- 
jectory at trap frequency 7^ = 2000, scattering length 
a — 0.1, and mean-field energy i?,„f = 2100. 

The excitation energy of the dynamically stabilized 
non-dissolving soliton in Fig 



11 



is A£' = E, 



mf 



111.0 which for a condensate with 10000 particles corre- 
sponds to a temperature of about T = 1.5 /iK. For the 
system with 7j, = 2000 and a = 0.1 stable oscillations 
can be found up to E^f = 2400. This is about 40 times 



the energy difference in Eq. (47) and yields for an esti- 
mation of the temperature T — 5.5 //K for N — 10 000 
particles. Calculations with a higher number of coupled 
Gaussians and with various initial conditions show that 
the oscillating soliton is stable. The choice of the ini- 
tial wave function is not crucial for the stability of the 
condensate. 
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FIG. 11. (Color online) Dynamics with 2 coupled Gaussians 
for trap frequency 7y = 2000, scattering length a = 0.1, and 
mean-field energy Emf ~ 2100. In (a) the parameters qx and 
Qz and in (b) the parameter Qy in the trap direction are shown 
as functions of time, (c) Projection of the trajectory in the 
(qxjQz) plane. The soliton does not dissolve although the 
mean-field energy is quite far above the threshold at = 



V. CONCLUSION AND OUTLOOK 

We investigated anisotropic quasi-2d solitons in dipo- 
lar Bose-Einstein condensates with the time-dependent 
variational principle using both the descriptive ansatz 
of a single Gaussian wave function and coupled Gaus- 
sians to calculate the ground state of the system. The 
lower limit of the energy obtained with a superposition 
of Gaussians is in full agreement with numerical grid cal- 
culations. For a given trap frequency we are able to 
give boundaries for the scattering length where stable 
solitons can exist. The energy gap between the mean- 
field energy of the ground state and the threshold energy 
where the soliton can dissolve is typically quite small. 
However, the investigation of the dynamics of the dipo- 
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lar BECs have revealed the existence of dynamicaUy sta- 
bihzed non-dissolving solitons at energies far above the 
threshold i?,„f ~ jy, and opens the possibility to cre- 
ate solitons at experimentally accessible temperatures [T] . 
This discovery may thus stimulate experiments on soli- 



tons in dipolar Bose-Einstein condensates. 
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